#pragma once

#include "local_search.h"

#include <limits>
#include <random>

namespace congruent_packing
{

template< class nesting_problem_t, template <class> class Objective >
class MSAS : public local_search_t
{
public:
	MSAS( nesting_problem_t & in_nesting_info )
	: nesting_info(in_nesting_info),
	  objective(in_nesting_info),
	  min_sigma(0.01),
	  Th(0.02),
	  annealing_parameter(0.98),
	  max_failure_count(3),
	  disp_factor(0.5)
	{}

	Solution work()
	{
		Solution & solution = nesting_info.points;

		Solution BSF = solution; // best solution found so far

		double sigmaI = 1.0;
		double min = -std::numeric_limits<double>::max();

		long int iteration = 0;

		size_t solutionSize = solution.size();

		const long int max_iterations = 1000000;

		while( sigmaI > min_sigma && iteration <= max_iterations )
		{
			for( size_t i = 0; i < solutionSize; i++ )
			{
				size_t failureCount = 0;
				size_t successCount = 0;

				double sigma = sigmaI;

				PointF point;

				while( failureCount < 4 * max_failure_count && successCount == 0 )
				{
					sigma = ( failureCount > max_failure_count ) ? disp_factor * sigma : sigma;

					point = nesting_info.points[ i ];

					point.x += gauss_distr(rand_engine,std::normal_distribution<>::param_type(0, sigma));
					point.y += gauss_distr(rand_engine,std::normal_distribution<>::param_type(0, sigma));

					nesting_info.container.clamp_point( point );
				
					if( objective( point, i, Th ) )
					{
						solution[i] = point;
						successCount++;
					}
					else
					{
						point = nesting_info.points[ i ];

						point.x -= gauss_distr(rand_engine,std::normal_distribution<>::param_type(0, sigma));
						point.y -= gauss_distr(rand_engine,std::normal_distribution<>::param_type(0, sigma));

						nesting_info.container.clamp_point( point );

						if( objective( point, i, Th ) )
						{
							solution[i] = point;
							successCount++;
						}
						else
						{
							failureCount++;
						}
					}
				}
			}

			double D = Mathematica::minDistance( solution );

			if( on_new_solution )
			{
				on_new_solution( iteration, D );
			}

			if( min < D )
			{
				min = D;

				if( on_best_record )
				{
					on_best_record( iteration, min );
				}

				BSF = solution;
			}

			iteration++;

			annealing_schedule( Th );
			annealing_schedule( sigmaI );
		}

		solution = BSF;

		return solution;
	}
	
	void annealing_schedule( double & temperature )
	{
		temperature *= annealing_parameter;
	}
		
public:
	double min_sigma;
	double Th;
private:
	double annealing_parameter;
	double disp_factor;
	size_t max_failure_count;

	std::normal_distribution<> gauss_distr;
	std::default_random_engine rand_engine;

	Objective<nesting_problem_t> objective;
	nesting_problem_t & nesting_info;
};

}